aov_duncan <- function(df, c1, c2) {

  model <- paste("y", paste(LETTERS[1:length(c1)], collapse = "*"), sep = "~")
  inter <- paste(LETTERS[1:length(c1)], collapse = ":")
  inter_name <- paste(names(df)[c1], collapse=":")

  res <- list()
  
  for(i in 1:length(c2)){
    df2 <- df[,c(c1,c2[i])]
    names(df2) <- c(LETTERS[1:length(c1)], "y")
    
    # user should adjust here
    fit <- aov(eval(parse(text = model)),data=df2)
    aov.res <- summary(fit) # list
    
    res1 <- list()
    fn <- names(df2)[1:length(c1)]
    for(j in 1:length(c1)){
      res1[[j]] <- duncan.test(fit,fn[j],alpha=.05)$groups
    }
    names(res1) <- fn
    
    res2 <- list()
    dfE <- aov.res[[1]]['Residuals','Df']
    MsE <- aov.res[[1]]['Residuals','Mean Sq']
    
    # user should adjust here
    res2[[inter_name]] <- with(df2, 
                               duncan.test(y,eval(parse(text = inter)),
                                           DFerror=dfE,MSerror=MsE))$groups
    
    res[[i]] <- list(aov.res,res1,res2)
    names(res[[i]]) <- c('aov','comp1','comp2')
  }
  names(res) <- names(df)[c2]
  return(res)
}